Multi-dimensional resilience: A quantitative exploration of disease outcomes and economic, political, and social resilience to the COVID-19 pandemic in six countries

The COVID-19 pandemic has highlighted a need for better understanding of countries’ vulnerability and resilience to not only pandemics but also disasters, climate change, and other systemic shocks. A comprehensive characterization of vulnerability can inform efforts to improve infrastructure and guide disaster response in the future. In this paper, we propose a data-driven framework for studying countries’ vulnerability and resilience to incident disasters across multiple dimensions of society. To illustrate this methodology, we leverage the rich data landscape surrounding the COVID-19 pandemic to characterize observed resilience for several countries (USA, Brazil, India, Sweden, New Zealand, and Israel) as measured by pandemic impacts across a variety of social, economic, and political domains. We also assess how observed responses and outcomes (i.e., resilience) of the COVID-19 pandemic are associated with pre-pandemic characteristics or vulnerabilities, including (1) prior risk for adverse pandemic outcomes due to population density and age and (2) the systems in place prior to the pandemic that may impact the ability to respond to the crisis, including health infrastructure and economic capacity. Our work demonstrates the importance of viewing vulnerability and resilience in a multi-dimensional way, where a country’s resources and outcomes related to vulnerability and resilience can differ dramatically across economic, political, and social domains. This work also highlights key gaps in our current understanding about vulnerability and resilience and a need for data-driven, context-specific assessments of disaster vulnerability in the future.


D
Associations between resilience metrics 13 E Additional modeling results 15

F Economic analysis 16
A Data sources 1 In this section, we provide additional information about the data sources and data 2 resolutions used. Table A.1 provides accession information for the data, including 3 website links. Table A.2 provides a catalog of the time and spatial resolution of each of 4 these variables in each of the countries. Although not shown here, some data sources 5 were available at further granular resolutions (e.g., daily cases or county-level data), but 6 all data sources were aggregated up to the weekly and administrative level 1 resolutions. 7 Resolutions in Table A  US unemployment data 08-03-2021 01/2010 to 06/2021 [11] World Input-Output Database (WIOD) Input-Output tables 6-16-2021 Pre-pandemic, 2016 release [12] 1 Time coverage corresponds to the longest time period with available data for any country. "Pre-pandemic" indicates a time-invariant variable collected for some pre-pandemic time. 1 Additional details about data sources can be found in Table A.1. 'Adm0' corresponds to (country-level) administrative level 0 resolution, and 'Adm1' corresponds to more granular (regional) administrative level 1 resolution. 2 Data are missing for at least one region.   In this section, we describe the strategy used to estimate COVID case under-reporting 13 in the main paper. We used a modification of the method proposed in Lau et al. 14 (2020) [13]. 15 Lau et al. (2020) [13] proposes estimating under-reporting based on the comparison 16 between the observed case-fatality rate (CFR) among reported COVID-19 cases and an 17 assumed "true" infection fatality rate (IFR) among all infected individuals. The first 18 challenge is in estimating the case fatality rate in the observed data, and the second 19 challenge is in defining the "true" infection fatality rate.

20
Several authors including Lau et al. (2020) [13] propose estimating the CFR as the total deaths by week w as a proportion of the total cases by week w-2. The goal of the time lag in deaths is to partially account for the delay in deaths among people diagnosed with COVID-19 in a given week. This strategy provides an estimated CFR over the entire course of the pandemic. However, since testing practices changed dramatically across the pandemic and our interest is in the under-reporting in a given time period rather than overall across the entire pandemic, we propose an alternative definition of the CFR within the last month as follows: CF Rrolling(w) = number of deaths between week w and w-4 number of cases between week w-2 and w-6 , (Eq. S1 ) where CF Rrolling(w) was then smoothed across 3 week intervals. The proposed rolling 21 monthly CFR estimates were designed to be used for correcting incident cases for 22 under-reporting rather than to correct cumulative cases, as described later on.

23
After we calculate the CFR, we need to define the target IFR that will serve as the 24 assumed "truth" for each administrative level 1 region. We will assume the true IFR by 25    We hypothesize that rate of vaccination in the population will have an impact on the 47 "target" IFR, even adjusting for age, and the method proposed in Eq. S2 will not 48 incorporate this. The complicated relationship between vaccination may be related to 49 the time since the last vaccination, receipt of a booster, the type of circulating COVID 50 variant, etc. Therefore, it is not trivial to fully account for the impacts of vaccination on 51 the target IFR.

52
As a rough sensitivity analysis, however, we propose a method for estimating under-reporting rates, assuming different simple relationships between vaccination and IFR. In particular, suppose vaccination impacts the likelihood of death among infected individuals multiplicatively, i.e., the IFR among vaccinated individuals is θ times the IFR among unvaccinated individuals. Here, we define vaccinated individuals as individuals reported to by fully vaccinated, with or without a booster. Then, we could calculate a rough vaccination-adjusted IFR estimate using the following equation: where p vax is the proportion of the population that is fully vaccinated, IFR refers to the 53 age-adjusted infection fatality ratio for unvaccinated individuals based on estimates 54 from O'Driscoll et al. (2021) [14], and θ ∈ (0, 1] is a sensitivity parameter controlling 55 the overall impact of vaccination on the IFR. Then, we obtain updated estimates of the 56 under-reporting factors by implementing the method in Eq. S2 using the new 57 vaccination-and age-adjusted estimate of the IFR.  66 We also notice that Sweden has estimated under-reporting factors less than 1 during 67 the spring/summer of 2021, as illustrated by Fig B.3. This indicates observed case 68 rates that are higher than expected, relative to the observed death rates. One possible 69 explanation for this is that COVID-related deaths may be under-reported. However, 70 this death under-reporting would need to be substantial to explain estimated 71 under-reporting factors of around 0.5. An alternative explanation is that the target IFR 72 values used in this analysis adjust for age and vaccination rates in the population on 73 average, but they do not adjust for other factors such as comorbidities, access to 74 healthcare, occupational exposures, etc. Failure to account for these additional 75 characteristics may result in poorly-specified "target" IFR values. This limitation may 76 make comparison of under-reporting rates between countries difficult; however, relative 77 changes in under-reporting factors within and between countries may be more reliable, 78 since many of these population characteristics related to IFR are time-invariant. In this section, we provide additional metrics describing the relationships between 90 resilience metrics over time. In Fig D.1, we provide the p-values related to pairwise 91 bivariate Granger tests as described in the main paper. The goal of these statistical 92 tests is to evaluate the relationship between a given predictor (columns) and 93 time-lagged past values of other predictors (rows). For all countries, past corrected case 94 counts are associated with future death rates. This is an artifact of the method used to 95 estimate the corrected cases, which directly involves use of death weeks two weeks in the 96 future. Fig D.       1 This table presents parameter estimates and 95% credible intervals from several varieties of mixed models for log-cases per 1000 for each of four countries separately. Model 1 contains a random intercept for region (administrative level 1) and an autoregressive error structure over time within each region. The column "Model 1a" provides results from fitting Model 1 to data where missing case counts (nearly always occurring prior to the first reported case for each region) are excluded. "Model 1b" provides results from fitting Model 1 to data where the missing cases count values are replaced with zeros and where government policy is excluded as a predictor from the model. "Model 2" uses the same data as "Model 1b" but also incorporates the spatial correlation between neighboring states within each country. Both Models 1 and 2 also include an 8-degree natural spline of calendar time. Population density (for Sweden) and governmental policy (all countries) were excluded from Model 2 due to missing data in these regions making the spatial correlation not able to be computed. In this section, we give a brief description of an Input-Output (IO) model. An IO model 101 is created starting from observations for a particular area, usually a nation or a region. 102 The economic activities of an area are classified into sectors/industries. We use the 103 International Standard Industrial Classification (ISIC) revision 3 industry classification 104 in this work, a standard defined by UN Statistical Division [15]. National data using 105 different classifications are mapped into ISIC code. See Table F.1 for the ISIC sectors' 106 definitions and Table F.2 for the ISIC sections, an aggregation of sectors.

107
The Input-Output relationships can be summarized as follows. Suppose the quantity x i is the total production of sector i and f i the total final demand of sector's i products.
In that case, we can represent the relationship between industries, i.e., how much one industry is selling to other sectors, by the following equation: where z i,j is the sale of sector i to sector j. We can express the relations between industries in matrix notation: x The total output of sector j is denoted x j . The ratio between z ij and x j is called the technical coefficient. The technical coefficient represents the amount of input from sector i expressed in dollars required to produce a one-dollar output in the industry j. We can rewrite F.1 as: If we group the x on the left side we have: x n = f n that in matrix notation corresponds to: For a given f , we solve for x: of COVID cases to determine the industry-specific COVID cases time series. 119 We use two data sources to compute the weights. The first one is published by the 120 California Department of Public Health [16]: following California rules, employers are 121 required to report workplace COVID-19 outbreaks to their local health department.

122
The data are published using the Census Industry Code.

123
The second data set is from the Employment Development Department of The State 124 of California, and it contains monthly employment data coded using the NAICS 125 industry classification. Using the "Census Industry Code List" published by Census [17] 126 we mapped the California outbreaks data to the NAICS industry classification. The two 127 datasets, outbreaks and employment, are finally mapped from NAICS to the ISIC 128 industry classification standard, the same used in the IO data [18]. From the reclassified 129 industry data we compute c i , the fraction of infected workers in industry i. The 130 industry infection weight r i is defined as: We assume that the estimated industry infection weight r i is the same across nations.

148
The sensitivity s i of the sector i to the workforce is the ratio between labor and the 149 total production requirements: s i = Li Li+Ki .

150
Using the national daily time-series data of the number of newly infected people, we 151 determine the number of infected people per industry sector i at time t, Ψ i (t) as: where δ is the fraction of working population over the total population, w i is the  s i , the sensitivity to the output of sector i to the labour factor: where p i is the production level of sector i. We assume that δ, p i , and s i are constant 159 during the analysis period. 160 Finally, we compute the economic drop in production d(t) for each sector, by propagating the shock across the economy of interconnected sectors through the input-output multipliers: where L is the Leontief Inverse matrix, cp(t) = (cp 1 (t), ..., cp n (t)) is the vector of the 161 change in production as in Eq. S6 . The drop in production, d(t), incorporates both the 162 direct and indirect effects of the shock on the economy.    1   2   3   4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38 39 40 41 424344 45 46 47 48  49  50  51  52  53  54   55   56   India   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38 39 40 41 424344 45 46 47 48  49  50  51  52  53  54 55 56 Figure F.3. Input-output tables network representation. The input-output table is pictured as a directed weighted graph. Each vertex represents an economic sector, and it is labeled with an index: see Table F.1 for a complete description of the sectors. The size of each vertex is proportional to the degree of the vertex, i.e., the number of its adjacent edges. The thickness of the graph edges is proportional to the connection strengths between sectors. The graph shows the most important connections between sectors: only the top 5 percent of the edges' weight distribution is visualized in the network.   Fig F.3), in the second column the economic sector identifier, and in the third column the sector definition.